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Abstract 

Presented  is  quantum  lattice-gas  model  for  simulating  the 
time-dependent  evolution  of  a  many-body  quantum  mechanical 
system  of  particles  governed  by  the  non-relativistic  Schroedinger 
wave  equation  with  an  external  scalar  potential.  A  variety  of 
computational  demonstrations  are  given  where  the  numerical 
predictions  are  compared  with  exact  analytical  solutions.  In 
all  cases,  the  model  results  accurately  agree  with  the  analytical 
predictions  and  we  show  that  the  model’s  error  is  second-order 
in  the  temporal  discretization  and  fourth-order  in  the  spatial 
discretization.  The  difficult  problem  of  simulating  a  system  of 
fermionic  particles  is  also  treated  and  a  general  computational 
formulation  of  this  problem  is  given.  For  pedagogical  purposes, 
the  two-particle  case  is  presented  and  the  numerical  dispersion 
of  the  simulated  wave  packets  is  compared  with  the  analytical 
solutions. 

KEY  WORDS:  Schroedinger  wave  equation,  quantum  comput¬ 
ing,  quantum  lattice  gas,  quantum  mechanics,  computational 
physics 


1  Introduction 

Feynman’s  work  regarding  quantum  mechanical  com¬ 
puters  used  to  simulate  physical  quantum  mechanical 
behavior  in  a  numerically  efficient  way  [1,  2,  3]  has  sub¬ 
sequently  led  to  several  series  of  research  papers  con¬ 
cerned  with  a  variety  of  details  involving  the  particu¬ 
lar  quantum  algorithm  that  best  represents  the  Feyn- 

*This  work  is  supported  by  the  Air  Force  Office  of  Scientific 
Research  under  the  Quantum  Computation  for  Physical  Mod¬ 
eling  initiative. 


man  path  integral  [4,  5].  The  quantum  algorithmic 
approach  is  based  on  a  two-component  complex  field 
defined  on  a  discrete  spacetime  lattice  where  unitary 
matrices  act  locally  on  the  field  causing  its  temporal 
evolution  in  discrete  time  steps.  Using  such  a  spatially 
discrete  field  makes  it  possible  to  computationally  rep¬ 
resent,  in  the  long  wavelength  limit  of  modes  in  the 
discrete  system,  the  dynamical  time-dependent  evolu¬ 
tion  of  a  continuous  wave  function  in  a  manner  that  is 
numerically  efficient. 

In  1994  Bialynicki-Birula  presented  a  general  quan¬ 
tum  algoithmic  approach  of  this  kind  for  modeling 
the  Weyl,  Dirac,  and  Maxwell  equations  on  a  body- 
centered  cubic  lattice  in  three-dimensions  [6].  In  a  se¬ 
ries  of  papers  on  simulating  the  one-dimensional  Dirac 
equation  [7,  8,  9],  Meyer  presented  a  quantum  algo¬ 
rithm  similar  to  that  of  Bialynicki-Birula  with  a  va¬ 
riety  of  numerical  simulations  including  the  effects  of 
boundary  conditions,  inhomogeneities,  and  an  exter¬ 
nal  scalar  potential.  Meyer  set  the  quantum  algorithm 
for  the  discretized  path  integral  in  the  context  of  what 
is  called  the  quantum  lattice  gas  method  and  his  algo¬ 
rithm  is  equivalent  to  the  one- dimensional  version  of 
the  Bialynicki-Birula  quantum  algorithm  for  the  Dirac 
equation. 

Contemporaneously  with  Meyer,  two  other  series 
of  papers  on  quantum  lattice-gas  models  of  the  one¬ 
dimensional  non-relativistic  Schroedinger  wave  equa- 
tiori  -'i^T-p  nrfisfinted  bv  Sued  and  Benzi  [10.  11]  and  by 
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Boghosian  and  Taylor  [12,  13,  14].  The  approach  taken 
by  Succi  and  Benzi  is  somewhat  more  computationally 
oriented  in  that  they  begin  with  a  “kinetic”  lattice 
Boltzmann  equation  of  motion1  (effectively  the  one¬ 
dimensional  Dirac  equation  in  the  Majorana  represen¬ 
tation)  and  show  that  the  Schroedinger  wave  equation 
emerges  as  the  governing  equation  of  motion  for  the 
slow  mode  in  the  long  wavelength  “hydrodynamic” 
limit.  That  is,  Succi  and  Benzi  observed  that  the 
“macroscopic  scale”  Schroedinger  wave  equation  arises 
from  the  “mesoscopic  scale”  Dirac  equation  in  a  man¬ 
ner  quite  analogous  to  how  the  macroscopic  Navier- 
Stokes  hydrodynamic  fluid  equation  arises  from  the 
mesoscopic  kinetic  Boltzmann  equation  through  the 
Chapman-Enskog  expansion. 

Boghosian  and  Taylor’s  approach  follows  more 
along  the  lines  of  Meyer’s  approach  in  that  their  model 
is  developed  as  a  generalization  of  the  classical  lattice 
gas  method.  A  kinetic  transport  equation,  now  for- 
muated  directly  at  the  “microscopic  scale,”  again  leds 
to  the  Schroedinger  wave  equation  in  the  continuum 
limit.  The  Boghosian  and  Taylor  quantum  lattice  gas 
model  focuses  on  solving  the  many-body  Schroedinger 
wave  equation  with  an  arbitrary  scalar  potential  in  an 
arbitrary  number  of  spatial  dimensions.  They  ana¬ 
lytically  argue  for  an  exponential  numerical  speedup 
arising  from  simulation  in  the  many-body  sector  of 
the  full  Hilbert  space  carried  out  simultaneously  us¬ 
ing  quantum  superposition  of  states.  The  Boghosian 
and  Taylor  version  of  the  quantum  algorithm  is  also 
cast  explicitly  for  direct  implementation  of  an  array  of 
quantum  bits  [13].  Polley  has  recently  presented  an  ar¬ 
gument  for  adding  both  an  external  scalar  and  vector 
potential  into  a  quantum  lattice-gas  model  by  analyt¬ 
ically  demonstrating  the  discrete  model’s  invariance 
with  respect  to  a  general  local  gauge  transformation 
[18]- 

A  characteristic  feature  of  all  these  quantum  algo¬ 
rithms,  used  to  model  the  dynamical  behavior  of  ei¬ 
ther  relativistic  or  non-relativistic  quantum  particles, 
is  that  the  governing  wave  function  is  well  approxi¬ 
mated  as  one  approaches  the  continuum  limit  where 
the  grid  resolution  of  the  spatial  lattice  become  in¬ 
finite  (the  lattice  cell  size  approaches  zero).  There¬ 
fore,  from  the  point-of-view  of  the  modeler,  there  ex¬ 
ists  a  “microscopic  scale”  where  the  unitary  dynamical 
rules  are  locally  applied  in  a  discretized  fashion  and 
time  advances  forward  in  incremental  units  in  a  way 
that  is  quite  artificial.  Yet  at  the  “macroscopic  scale,” 
which  corresponds  to  the  long  wavelength  limit  of  the 

1  The  classical  lattice  Boltzmann  equation  was  popularized 
with  its  application  to  computational  fluid  dynamics  [15, 16, 17]. 


dynamical  modes  in  the  discrete  system,  there  is  an 
emergent,  effective  field  theory  for  a  complex  ampli¬ 
tude  field,  continuous  and  differentiable  in  both  space 
and  time,  that  exactly  obeys  the  physical  quantum  me¬ 
chanical  equations  of  motion.  At  the  small  scale  (char¬ 
acterized  by  the  lattice  cell  size)  one  imagines  a  sys¬ 
tem  of  fermionic  particles  undergoing  local  collisions 
and  translation  to  nearby  nodes  of  the  lattice.  Each 
of  these  fermionic  particles  occupies  a  local  positional 
state  at  a  specific  lattice  node  with  a  certain  probabil¬ 
ity  amplitude.  All  the  possible  locations  of  the  actual 
physical  quantum  particle  are  effectively  modeled  by 
the  interfering  set  of  probability  amplitudes  associated 
with  this  system  of  fermionic  particles.  That  is,  all  the 
possible  pathways  are  modeled  simultaneously  using  a 
kind  of  kinetic  system  of  locally  interacting  fermions 
on  the  small  scale. 

Seen  as  a  kinetic  system  then,  we  may  expect  that 
there  exists  a  local  equilibrium  configuration  of  parti¬ 
cles.  We  require  that  this  configuration  be  an  eigen- 
ket,  with  unity  eigenvalue,  of  the  local  unitary  colli¬ 
sion  operator  that  is  uniformly  and  spatially  homo¬ 
geneously  applied  to  the  lattice-based  two-component 
field.  Then  the  dynamical  lattice-based  system  on  all 
lattice  nodes  undergoes  local  relaxation  towards  this 
equilibrium  configuration.  However,  unlike  a  classi¬ 
cal  kinetic  system,  the  global  configuration  of  particles 
does  not  relax  to  a  single  steady-state  equilibrium.  In¬ 
stead,  there  are  many  steady-state  global  equilibrium 
configurations  which  effectively  are  the  energy  eigen¬ 
states  of  the  “macroscopic  scale”  quantum  mechanical 
equation  of  motion  modeled  by  the  quantum  lattice 
gas.  If  the  quantum  lattice  gas  system  is  initialized  in 
any  one  of  the  energy  eigenstate  global  configurations, 
the  macroscopic  scale  configuration  of  the  system  will 
remain  fixed  in  time  albeit  the  microscopic  configu¬ 
ration  of  particles  would  be  continually  changing  at 
every  unit  time  step.2  In  the  end,  the  Feynman  path 
integral  is  efficiently  and  accurately  recast  as  a  kinetic 
dynamical  process  computed  in  parallel  efficiently  on 
a  spacetime  lattice. 

In  this  paper  we  do  not  argure  that  the  quantum 
lattice  gas  dynamical  rules  represent  a  local  hidden 
variable  theory  of  quantum  mechanics.  Although  fun¬ 
damental  arguments  can  be  made  to  limit  the  possi¬ 
ble  form  of  the  local  unitary  operator  in  the  quan¬ 
tum  lattice  gas  [6,  7],  these  arguments  lead  to  an  algo¬ 
rithm  suited  for  implementation  on  a  quantum  com- 

2Given  a  finite  size'lattice  used  for  modeling  purposes,  each 
local  configuration  oscillates  in  time  even  when  the  global  con¬ 
figuration  is  a  time-independent  energy  eigenvalue.  However, 
the  amplitude  of  the  oscillation  does  approach  zero  as  the  lat¬ 
tice  size  becomes  infinite. 
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puter.  In  our  investigation  of  a  suitable  local  unitary 
collision  operator  we  have  found  that  one  quantum 
gate  in  particular,  the  square-root-of-swap  given  be¬ 
low  in  Section  2.3,  is  especially  useful  for  modeling 
the  Schroedinger  wave  equation  in  that  the  local  equi¬ 
librium  configuration  discussed  above  is  an  eigenket 
of  this  gate  and  with  unity  eigenvalue.  Therefore,  we 
have  selected  the  square-root-of-swap  gate  as  our  basic 
model  quantum  gate.  As  demonstrated  in  Section  2.6, 
we  find  that  this  quantum  gate  leads  to  an  overall  mod¬ 
eling  error  that  is  second-order  in  the  temporal  dis¬ 
cretization  and  fourth-order  in  the  spatial  discretiza¬ 
tion.  Another  point  regarding  this  particular  gate  is 
that  when  measurements  are  periodically  made  of  the 
state  of  qubits  in  the  system,  which  destroys  quantum 
superpositions  and  entanglements  in  the  system,  the 
macroscopic  scale  behavior  of  the  quantum  lattice  gas 
system  is  governed  by  the  classical  diffusion  equation 

[19]- 

This  paper  is  organized  as  follows.  First,  in  Sec¬ 
tion  2  we  describe  the  basic  quantum  algorithm,  in 
particular,  how  we  encode  the  wave  function,  how  we 
use  two  basic  quantum  gates  applied  to  the  qubits  at 
the  nodes  of  the  lattice.  We  first  describe  the  algo¬ 
rithm  using  matrices  and  then  we  describe  an  equiv¬ 
alent  finite  difference  formulation.  Using  the  micro¬ 
scopic  finite  difference  equations,  we  then  derive  the 
effective  field  theory  at  the  macroscopic  scale  where 
both  the  lattice  cell  size  and  the  update  time  step  ap¬ 
proach  zero.  Diffusive  ordering  holds,  as  is  typical  for 
lattice  gas  systems,  where  small-scale  temporal  fluc¬ 
tuations  in  the  wave  function  go  as  the  square  of  the 
magnitude  of  the  small-scale  spatial  fluctuations.  To 
confirm  that  our  derivation  is  correct  and  to  test  the 
validity  of  the  quantum  lattice  gas  model,  we  test  the 
time-dependent  dispersion  of  a  free  Gaussian  packet. 
We  also  test  the  system  when  it  is  initialized  in  an  en¬ 
ergy  state,  which  is  a  fixed  macroscopic  scale  steady- 
state  configuration.  We  find  that  as  we  halve  the  lat¬ 
tice  cell  size  (double  the  spatial  resolution),  the  cumu¬ 
lative  error  in  the  model  drops  by  a  factor  25  45  =  43.7. 

Second,  in  Section  3  we  show  how  an  external  scalar 
potential  may  be  modeled  in  the  quantum  lattice  gas 
system  by  inducing  a  local  phase  rotation  in  the  qubits 
at  each  node  of  the  lattice.  The  qubits  at  a  lattice  node 
are  phase  rotated  by  an  amount  corresponding  to  the 
strength  of  the  spatially  dependent  external  potential 
at  that  node.  We  show  how  this  local  phase  rotation,  a 
kind  of  gauge  transformation,  produces  an  additional 
potential  energy  term  in  the  Schroedinger  wave  equa¬ 
tion.  We  then  test  the  quantum  algorithm  against  two 
well-known  cases  of  harmonic  oscillation  in  a  parabolic 


potential  well  and  quantum  scattering  off  of  and  tun¬ 
neling  through  a  constant  potential  energy  barrier.  In 
both  cases,  the  model  behaves  as  expected. 

Third,  in  Section  4  we  test  how  well  the  quantum 
lattice  gas  can  simulate  the  simultaneous  dispersion 
of  two  fermionic  particles.  In  the  case  of  multiple 
particles,  the  operational  quantum  gate  sequence  to 
handle  the  many-body  case  is  identical  to  the  sin¬ 
gle  particle  case  presented  in  Section  2.  Therefore, 
the  implementation  of  either  situation  on  a  quan¬ 
tum  computer  would  be  identical  as  well,  except  for 
state-preparation  and  measurements.  However,  since 
at  present  no  quantum  computer  exists  that  can  test 
the  algorithm  presented  here,  we  are  forced  to  con¬ 
sider  the  implementation  on  a  classical  computer  and 
here  the  implementation  for  a  many-body  problem  is 
much  more  complex  than  for  the  single-body  prob¬ 
lem.  Nevertheless,  we  present  a  general  formulation 
of  the  quantum  gates  in  a  second  quantized  represen¬ 
tation  where  the  basic  computational  operations  are 
creation  and  annihilation  of  local  particle  occupancies. 
The  advantage  of  such  an  implementation  is  that  it 
is  straightforward  to  implement  the  fundamental  cre¬ 
ation  and  annihilation  operation  is  a  way  that  respects 
the  anti-commutation  relations  for  fermionic  particles. 
We  demonstrate  that  the  macroscopic  scale  behavior 
of  the  quantum  lattice  gas  agrees  with  the  exact  time- 
dependent  solution  of  the  two-body  wave  equation. 

Finally,  we  conclude  this  paper  with  a  short  de¬ 
scription  of  results  and  share  some  lessons  we  learned 
after  using  the  quantum  lattice  gas  model  extensively. 
We  also  point  out  some  future  directions  that  may  be 
taken  to  expand  the  usefulness  of  the  model. 

2  Quantum  Algorithm  for  a  Sin¬ 
gle  Free  Particle 

We  describe  the  quantum  lattice-gas  algorithm  for 
modeling  the  Schroedinger  wave  equation  by  consid¬ 
ering  the  simplest  case  of  a  single  free  particle  in  a 
one-dimensional  space.  In  this  simple  case,  the  wave 
function  ip(x,  t )  obeys  the  following  partial  differential 
equation  in  the  position  representation 

,tdip(x,t)  _  h2  d2ip(x,t) 

dt  2m  dx 2  ’  W 

where  h  is  Planck’s  constant  and  m  is  the  mass  of 
the  quantum  particle.  Here  ip(x,  t)  is  a  continuous 
probability  amplitude  field  (e.g.  a  continuous  complex 
field). 
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2.1  Encoding  the  Wave  Function 

To  “program”  a  quantum  computer  to  simulate  (1), 
it  is  necessary  to  first  formulate  an  encoding  scheme 
where  a  collection  of  qubits  is  used  to  store  the  value 
of  the  wave  function.  Since  the  number  of  qubits  in 
any  quantum  computer  is  necessarily  a  finite  number, 
the  wave  function  will  have  to  be  approximated  in  the 
usual  way  by  discretizing  a  physically  continuous  am¬ 
plitude  field  into  an  artificially  discrete  and  finite  set  of 
complex  numbers.  To  do  this,  let  us  begin  with  a  one¬ 
dimensional  spatial  lattice  with  L  number  of  nodes. 
With  each  node  of  the  lattice  we  associate  a  position 
basis  ket  denoted  by  |*j),  where  0  >  l  >  L  —  1.  The 
discretized  system  ket  in  the  position  basis  is 

L—l 

W)  =  (2) 

1=0 

where  c;  =  (aq|i/>)  is  a  complex  number.  In  other 
words,  the  basic  approach  to  model  the  single  par¬ 
ticle  wave  function  governed  by  (1)  is  to  express  \ip) 
as  a  stun  of  all  the  possible  ways  the  particle  can  be 
situated  on  the  lattice  with  a  probability  amplitude  c; 
associated  with  each  possible  location  \xi). 

In  our  model,  we  assign  two  qubits  to  each  node  of 
the  lattice,  for  a  total  of  2 L  qubits  in  the  whole  quan¬ 
tum  computer.  The  qubits  that  reside  at  the  Zth  node 
of  the  lattice  are  denoted  by  |§q)  and  \q[)  and  they  are 
used  to  encode  the  coefficient  c;  of  (2)  of  the  position 
ket  for  that  node.  Each  qubit  is  a  two-level  quantum 
system  |q£)  =  ai|0)+/?£|l),  where  \ala\2 +  \/3la\2  =  1  for 
a  =  0  or  1  and  0  >  l  >  L  —  1.  We  consider  each  qubit 
to  be  a  container  that  may  or  may  not  be  occupied  by 
the  quantum  particle.  The  quantum  particle  is  said  to 
occupy  the  ath  local  state  at  position  xi  when  /?*  =  1. 
Similarly,  the  ath  local  state  at  position  xi  is  said  to 
be  empty  when  /?*  =  0. 

To  see  how  the  qubit  encoding  works,  we  write  | ip) 
in  the  number  representation.  In  the  number  repre¬ 
sentation,  each  basis  state  is  expressible  as  the  ket 
lnoninonlnoni ' '  ‘noni),  where  nla  —  0  or  1  for  all  l 
and  a.  The  Boolean  variables  nla  are  called  the  number 
variables  and  they  correspond  to  a  binary  indexing  of 
the  basis  states  in  the  number  representation.  Since 
we  are  concerned  with  modeling  the  one-particle  wave 
equation,  we  need  consider  only  a  subset  of  all  the  ba¬ 
sis  states  where  only  one  of  the  number  variables  is 
1  and  all  the  rest  are  0.  This  subset  of  all  the  ba¬ 
sis  states  is  called  the  one-particle  sector.  There  are 
2 L  such  combinations  and  we  shall  label  these  with 
the  binary  encoding  formula  |22i+°),  for  a  =  0, 1  and 
0  >  l  >  L— 1.  Therefore,  the  system  ket  in  the  number 


representation  can  be  written  as 

i^  =  EE^+“i22i+°)>  ‘  .  (3) 

l— 0  o=0 

where  each  £2 i+a  is  a  probability  amplitude  (e.g.  com¬ 
plex  number). 

Now  for  each  position  ket  |aq)  there  are  two  cor¬ 
responding  basis  states  in  the  number  representation 
|22*)  and  |22I+1).  There  are  two  interfering  possibilities 
for  a  particle  to  occupy  the  1th  position  on  the  lattice. 
Therefore,  the  occupancy  probability  of  the  1th  node  is 
computed  by  first  summing  the  probability  amplitudes 
of  these  corresponding  basis  states  and  then  comput¬ 
ing  the  square  of  the  absolute  value  thereof.  In  other 
words,  the  coefficient  c;  in  (2)  is  set  equal  to  the  sum 
of  the  on-site  coefficients  in  (3) 

ci  =  £21  +  fri+i-  (4) 

The  definition  (4)  is  an  essential  part  of  the  quantum 
lattice-gas  model  presented  in  this  paper.  In  the  sec¬ 
tion  below,  where  we  analytically  predict  an  effective 
field  theory  for  our  artificially  discretized  model,  we 
explain  why  we  need  to  make  this  assignment.  We  will 
find  that  (4)  is  needed  for  the  predicted  effective  field 
theory  to  accurately  approximate  the  Schroedinger 
wave  equation  in  the  long- wavelength  limit,  which  is 
also  defined  below. 

2.2  Formulating  a  Suitable  Gate  Se¬ 
quence 

We  shall  require  that  the  algorithmic  scheme  be  at 
least  second  order  convergent  in  space,  so  that  as  we 
double  the  grid  resolution  (e.g.  double  the  number 
of  qubits)  we  in  turn  reduce  the  numerical  error  due 
to  the  field  discretization  by  a  factor  of  one-quarter. 
With  this  type  of  convergence  characteristic,  we  are 
assured  that  we  can  simulate  a  wave  function  gov¬ 
erned  by  the  Schroedinger  wave  equation  (1)  to  any 
arbitrary  degree  of  accuracy.  After  we  formulate  our 
algorithmic  scheme,  we  will  then  a  posteriori  verify  by 
direct  numerical  simulation  that  it  is  indeed  at  least 
a  second-order  convergent  numerical  scheme.  In  fact 
in  Section  2.6  we  will  find  that  our  numerical  scheme 
is  foUrth-order  convergent  with  an  error  that  goes  as 
(5x)4. 

To  simulate  the  quantum  behavior  of  the  wave  func¬ 
tion,  we  seek  to  develop  a  sequence  of  2-qubit  gate  op¬ 
erations  that  will  act  on  a  large  collection  of  qubits  in 
the  simplest  way.  We  impose  the  following  four  sim¬ 
plifying  constraints: 
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1.  All  quantum  gate  operations  are  homogeneous 
and  independent  of  space  and  time.  • 

2.  Only  a  single  quantum  gate  is  used  to  evolve  the 
wave  function  and  this  gate  is  applied  to  each  lat¬ 
tice  node  independently  (locality). 

3.  To  provide  communication  channels  between  lat¬ 
tice  nodes,  only  the  simplest  gate  is  used  (e.g.  a 
swap  gate). 

4.  Because  the  final  value  of  the  computed  wave 
function  depends  on  summing  interfering  possibil¬ 
ities  according  to  (4),  we  shall  use  the  Haddamard 
gate  at  the  very  end  of  the  simulation  prior  to 
making  a  measurement  of  the  wave  function  so 
that  a  single  qubit  at  each  node  will  encode  the 
probability  amplitude  ci  =  (xi\tp)  in  (2). 


The  full  collision  operator,  denoted  C,  which  acts 
on  the  system  ket  |?/’}  is  formed  by  a  L-fold  tensor 
product  over  the  local  collision  operators  U  applied 
homogenously  and  independently  on  each  node  of  the 
lattice 

L-l 

C=®U.  (7) 

1=0 

Let  us  denote  the  swap  operator  by  x im,v,  where 
H  and  v  index  any  two  qubits  in  the  system.  The 
streaming  operator,  denoted  Si,  causes  a  global  shift 
to  the  right  of  the  first  qubit  on  all  the  lattice  nodes. 
Therefore,  Si  can  be  represented  by  a  sequence  of  swap 
operators  acting  on  nearest  neighbors 

L—l 

2 

Si  =  U  X2l,2l+2-  (8) 

1=0 


With  two  qubits  per  node,  there  are  four  on-site 
basis  kets,  |0)  ®  |0)  =  (1,0, 0,0),  |0)  ®  |1)  =  (0, 1, 0,0), 
|1)®|0>  =  (0,0, 1,0),  and  |1>  ®  |1>  =  (0,0, 0,1).  In  the 
context  of  a  quantum  lattice-gas  model,  the  unitary 
matrix  U  is  called  the  local  collision  operator  and  the 
on-site  ket  \u)  =  |0)®|1)-|-|1)®|0)  =  (0,1, 1,0)  is  called 
the  number  density  ket.  To  have  a  well  behaved  local 
equilibrium  associated  with  the  collision  process,  the 
local  collision  operator  must  have  the  number  density 
ket  as  an  eigenvector  with  unity  eigenvalue. 


2.3  Matrix  Representation 

The  quantum  gate  that  we  use  to  evolve  the  wave  func¬ 
tion,  which  is  applied  independently  on  a  site-by-site 
basis,  is  the  square-root-of-swap  gate 


The  reason  for  calling  this  the  square-root-of-swap 
gate  is  that  U 2  is  the  swap  gate  itself 


U-U  = 


(1  0  0  0\ 
0  0  10 
0  1  0  0  I 
0  0  0  1/ 


(6) 


The  two  nontrivial  eigenvalues  of  U  are  Ai  =  1  and 
A2  =  — t,  with  eigenvectors  |iq)  =  (0, 1, 1, 0)  and  ^2)  = 
(0,  —1, 1, 0),  respectively.  Also,  since  (5)  causes  mixing 
only  between  the  single-particle  basis  kets  |0)®|1)  and 
jl)  ®  |0),  it  conserves  particle  number.  So  (5)  is  an 
appropriate  choice  for  the  local  collision  operator. 


In  matrix  form,  x  is  a  22  x  22  permutation  matrix 
1  0  0  0  \ 


000  -1/ 

The  algorithm  we  use  to  model  the  Schroedinger  wave 
equation  involves  multiple  applications  of  the  collision 
operator  interleaved  with  streaming  operations  as  fol¬ 
lows:  5 

m+T/2))=Ei\m),  (10) 

where  the  square  root  of  the  evolution  operator  is 

e}  =  SjcSiC.  (11) 

Here  Sf  denotes  the  transpose  of  S±  and  is  the  inverse 
of  Si.  Application  of  Sj  causes  a  global  shift  to  the 
left  of  the  first  qubit  on  all  the  lattice  nodes.  One  full 
time  step  of  the  evolution  is 

\iJ(t  +  r))  =  Ei\ij)(t)).  (12) 

We  use  four  applications  of  the  collision  operator  in 
E\  because  C14  is  the  identity  operation.  Note  that 
Si  and  C  do  not  commute,  otherwise  (12)  would  be  a 
trivial  evolution  equation. 

Note  that  in  (12),  our  choice  of  streaming  the  first 
qubit  was  arbitrary.  A  more  balanced  algorithmic  ap¬ 
proach  would  treat  both  qubits  identically.  Therefore, 
we  could  alternatively  define  one  full  time  step  as 

m  +  T))  =  E2Ei\iP(t)),  (13) 

where 

El  =  SICS2C,  (14) 


0  0  10 
0  10  0 
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and  where  the  streaming  operator  S-2  causes  a  global 
shift  to  the  right  of  the  second  qubit  on  all  the  lat¬ 
tice  nodes.  The  advantage  of  using  the  balanced  al¬ 
gorithm  (13)  is  that  its  error  is  fourth-order  in  space 
whereas  for  the  unbalanced  algorithm  (12)  it  is  only 
third-order. 

2.4  Finite  Difference  Formulation 

It  is  possible  to  specify  the  quantum  algorithm  to 
model  the  Schroedinger  equation  without  the  use  of 
matrices.  Instead  we  can  write  down  a  set  of  finite 
difference  equations,  which  are  equivalent  to  (10),  but 
perhaps  simpler  to  comprehend  at  first  glance.  To  do 
this,  let  us  introduce  a  new  notation  for  the  2 L  prob¬ 
abilities  amplitudes  £2 i+a  in  (3).  We  will  denote  the 
two  complex  numbers  per  lattice  node  by  <po(xi,  tn) 
and  ipi(xi,tn).  That  is,  we  have  L-pairs  of  complex 
numbers.  Then,  the  quantum  algorithmic  operations 
(14)  can  be  expressed  as  follows: 


if  mod(n,  4)  =  0 

(15) 

^n) 

=  A*if0(xl,tn-1)  +  Atpi(xi,tn-i) 

(plipU  tn) 

=  Aipo{xi,tn-i)  +  A*tpi(xutn-i), 

if  mod(n,  4)  =  1 

(16) 

<Po(xU  tn) 

(po(xi-i,tn-i) 

<Pl  (p'h  tn) 

if  mod(n,  4)  —  2 

(17) 

tn) 

=  A*ipo(xi,tn-i)  +  A<pi{xi,tn-i) 

(plip^U  tn) 

—  Aipo(xi,tn-i)  +  A  <f\{xi, tn~ i), 

and  if  mod(n,4)  =  3 

(18) 

^PoiptU  tn) 

=  9^0(^1+11^71—1) 

<Pl  0 tn) 

(fi(xi,  tn— 1), 

where  A  =  \Jr\-  The  finite-difference  equation  pair 

(15)  is  equivalent  to  the  local  collision  operation  C, 
as  is  the  pair  (17).  The  equation  pairs  (16)  and  (18) 
are  equivalent  to  the  streaming  operations  S  and  ST, 
respectively.3 

This  finite-difference  representation  of  the  al¬ 
gorithm  is  nearly  identical  to  that  presented  by 
Boghosian  and  Taylor  in  1997  [14]  where  the  two  on¬ 
site  qubits  are  simultaneously  streamed  to  the  left  and 

3  Noting  that  A  +  A*  =  1,  this  set  of  finite  difference  equa¬ 
tions  can  be  expressed  in  a  more  compact  way 

VoOq+oWl)  =  Vo{xi,tn)  +  fio 

fli^lAn+l)  =  <pi(K;,tn) +fli 

where  e  =  (— l)n  and  fio  =  A(ipi  —  ipo)  and  fix  =  —Ho,  which 
has  the  standard  form  of  a  lattice-gas  transport  equation. 


right  after  collision  operation 

<Po(xl+1,tn)  =  A*ip0(xi,tn_1)  -  A<pi{xi,tn..i) 
<Pl(xi-i,tn)  =  -A<po(xi,tn-l)  +  A*(pi(xi,tn-1). 

■  (19) 

They  noted  that  after  four  time  steps,  the  total  am¬ 
plitude  ip(xi,tn)  =  <f>o(xi,tn)  +  4>i(xi,tn)  satisfies 
a  finite-difference  equation  which  approximates  the 
Schroedinger  equation  in  the  continuum  limit.  The 
two  essential  differences  between  the  improved  algo¬ 
rithm  (15)  through  (18)  presented  in  this  paper  and 
the  quantum  algorithm  (19)  appearing  in  [20]  is  that 
we  have  alleviated  the  problem  of  the  occurrence  of 
two  non-interpenetrating  lattice-gas  systems  indepen¬ 
dently  evolving  on  different  checker-board  sub-lattices 
and  we  have  doubled  the  numerical  accuracy.  This 
is  a  problem  that  occurs  when  both  on-site  qubits 
are  simultaneously  streamed  because  streaming  only 
a  single  qubit  at  a  time,  as  was  done  for  the  quantum 
lattice-gas  model  of  the  diffusion  equation  [19],  causes 
interactions  between  all  the  qubits  at  each  time  step. 

2.5  The  Governing  Partial  Differential 
Equation 

It  is  straightforward  using  a  symbolic  mathematics 
program,  and  tedious  by  hand,  to  use  the  update  rules 
(15)  through  (18)  to  algebraically  determine  the  value 
of  <po  and  <p\  at  a  later  time.  With  the  initial  wave 
function  set  at  to,  one  complete  cycle  of  the  algorithm 
is  completed  at  ts  (that  is,  tg—to  =  <5t).  With  the  wave 
function  defined  as  ip(xi,tn)  =  tpo{xi,tn)  +  <fii(xi,tn), 
the  result  after  one  cycle  is4 

1  -\.  i 

^{xi,ts)  = - —  ip(xi,to)  +i>(xi+1,to)  +  r/)(xi-i,to) 

[V'(a;(+2,  to)  +  ip(xi-2,  to)] ,  (20) 

Note  that  (20)  is  the  simplified  form  of  the  finite- 
difference  equation  at  the  macroscopic  scale  when  the 
system  is  very  close  to  local  equilibrium  throughout 
the  course  of  the  evolution  as  ipo(x,t)  =  ~ 

\'<p(x,  t )  for  all  x.  The  full  finite-difference  equation  is 
too  long  to  present  here,  but  is  given  in  Appendix  A. 
This  result  is  a  finite-difference  equation  for  the  fol- 

4  Note  that  the  result  (20)  is  accurate  up  to  fourth  order 
in  5x  only  in  the  situation  where  the  initial  system  is  in  local 
equilibrium  defined  by  (po(^l,tn)  =  In  the  more 

general  situation  when  the  system  is  not  in  local  equilibrium 
where  ipo(%l>tn)  A  <Pl(xi,tn)i  the  result  (20)  is  accurate  only 
up  to  third  order  in  Sx. 
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lowing  partial  differential  equation  governing  the  con¬ 
tinuous  amplitude  field  ?/’(a:,  t ) 


dip(x,t) 

dt 


+  0(8t 2 


1  Sx 2  d2ip(x,t) 

2  5t  dx 2 


+  0(5x 4),  (21) 


which  is  an  approximation  of  (1)  where  the  diffusion 
constant  is  h/m  =  5x2/5t  and  where  5x  is  the  lattice 
cell  size. 

If  one  adds  a  phase  angle  £  to  the  off-diagonal  com¬ 
ponents  of  collision  operator  (5)  to  obtain  a  slightly 
more  general  collision  operator 


/I 

0 

0 

0  > 

0 

1  i 

2  2 

(l  +  l)eiC 

0 

0 

(1  +  1)  ^iC 

1  i 

2  2 

0 

^0 

0 

0 

-1 ) 

then  the  resulting  governing  partial  differential  equa¬ 
tion  will  have  its  transport  coefficient  dependent  on 
this  phase  angle  as  follows: 


dip(x,t) 

dt 


+  0(St2)  = 


i  8x2  d2ip(x,t) 
2  sec  C  8r  dx 2 


+  0{SxA). 


(23) 

This  allows  us  to  simulate  a  quantum  system  where 
a  particle’s  mass  can  be  arbitrarily  large  m  —  sec  (, 
but  has  a  minimum  of  one.  Note  that  in  this  case  the 
error  is  cubic  and  is  proportional  to  sin  (.  So  for  very 
large  masses,  the  accuracy  of  the  model  is  reduced 
to  third-order  in  space.  Note  that  (23)  is  is  valid  ef¬ 
fectively  field  theory  at  the  macroscopic  scale  when 
the  system  is  very  close  to  local  equilibrium  where 
<Po(x,  t)  =  ipi{x,t)  ~  \^{x,t)  for  all  x. 


2.6  Numerical  Confirmations 

To  numerically  test  that  the  quantum  algorithm  (12) 
is  indeed  equivalent  to  the  finite-difference  equation 
(20)  and  to  see  just  how  good  of  an  approximation  of 
the  single-particle  Schroedinger  equation  it  is,  we  have 
preformed  two  simulations. 

In  the  first  simulation,  we  test  the  numerical  time 
evolution  of  a  Gaussian  packet 

=  —t L—  e~S^,  (24) 

a  2  7r2 

where  l  >  x  >  L  for  a  lattice  of  size  L  =  64£  and  where 
the  packet  width  is  a  —  L/ 10  as  shown  in  Figure  1. 

The  exact  anaytical  solution  of  (21)  is  obtained  by 
computing  the  Fourier  components  of  the  energy  basis  ' 
functions 

i  /-b/2 

a0  =  —  /  ij)(x,0)dx  (25) 

L  J-L/2 


10  20  30  40  50  60 

Position 


Figure  1:  Time  evolution  of  a  Gaussian  packet  for  a  single 
quantum  particle  overplotted  in  succession  where  the  x-axis  is 
the  position  on  a  64-node  lattice  in  units  of  the  lattice  spacing 
t  and  the  y-axis  is  the  probability  density  \ip(x,t)\2.  The  solid 
curves  are  the  exact  analytical  solution  and  the  circles  are  the 
data  from  the  quantum  lattice-gas  simulation  (the  initial  wave 
function  was  normalized,  therefore  the  area  under  each  curve  is 
one).  The  lattice  size  is  L  =  64f.  The  initial  Gaussian  packet  of 
with  a  =  L/ 10  at  t  =  0  is  centered  at  x  =  32 i  and  the  disper¬ 
sion  is  evident  by  observing  the  wave  function  at  the  subsequent 
times  t  =  50 r,  100r,  150r,  and  200r.  Periodic  boundary  condi¬ 
tions  were  used  and  nmBx  =  20  energy  eigenmodes  were  used 
to  generate  the  exact  solutions.  A  time  scale  factor  ts  =  1.04 
was  used  to  improve  the  agreement  between  the  numerical  and 
analytical  solutions. 


2 

/•i/2 

=  T 
A  j 

/  ip(x,  0)cos 

'-L/2 

(2n  *l) 

With  h  =  1  and  m=  1,  the  energy  eigenvalues  are 


2n2ir2dx2 
=  L26t 


(27) 


and  the  time-dependent  solution  to  (21)  plotted  in  Fig¬ 
ure  1  is 

^exactOM)  =  do  +  an  cos  (2mr— J  e~lEnt/t‘ .  (28) 

n=l 

Note  that  in  (28),  time  is  scaled  by  a  factor  ts  to  ac¬ 
count  for  kinetic  corrections  to  the  time  step.  As  the 
number  of  lattice  nodes  becomes  large,  this  scaling 
factor  approaches  one. 

The  second  test  of  the  quantum  lattice-gas  algo¬ 
rithm  as  a  model  of  the  Schroedinger  wave  equation  is 
the  measurement  of  its  numerical  convergence.  Mul¬ 
tiple  simulations  (10  in  total)  were  carried  out  for 
lattice  sizes  ranging  from  L  =  8£,  16£,  32£, ...  up  to 
L  =  8192f.  In  each  case  the  initial  state  of  the  simula¬ 
tion  was  the  ground  state  (a  sinusoidal  energy  eigen- 
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Figure  2:  Log-log  plot  of  the  numerical  error  versus  re¬ 
solving  grid  cell  size,  Sx ,  indicating  the  convergence  property 
of  the  quantum  lattice-gas  algorithm  (12)  and  (13)  for  the 
Schroedinger  equation.  The  data  (black  circles)  are  taken  from 
numerical  simulations  with  grid  sizes  from  L  =  8^  up  to  8192£ 
after  a  single  time  step  T  =  r.  The  solid  curves  are  best-fit 
linear  regression  with  a  slope  of  3.48  and  5.45  for  the  models 
defined  by  (12)  and  (13),  respectively.  These  results  demon¬ 
strate  third-order  and  fourth-order  convergence  in  space  for  the 
two  models,  respectively. 


state) 


tp(x,  0)  =  tp‘XBCt(x)  = 


cos  (2-kx/L) 

VlJ 2 


(29) 


Each  simulation  was  run  for  one  time  step  T  —  r  and 
the  numerical  error,  denoted  e,  from  the  exact  solution 
was  then  measured  using  the  following  formula 


e(L) 


1 


'\ 


£{|^T)|2-|^(*)|2}2.  (30) 


X=1 


We  define  the  grid  resolution  as  the  inverse  of  the  total 
number  of  lattice  points.  That  is,  for  a  box  of  size  1, 
the  resolving  cell  size  is  defined  as  5x  =  ■£.  A  plot 
of  the  error  versus  the  resolution  is  given  in  Figure  2. 
As  the  resolution  is  increased,  the  error  drops  off  as 
e(L)  ~  L~5A5. 


3  Adding  an  External  Scalar 
Potential 


e-iV^^(xl+1,to) 
1  —  i 


+  e 


[e-iV^)st^(Xl+2,to) 
o)]  . 


If  we  expand  the  potential  terms  in  the  arguments  of 
the  exponentials 


V(xi+i)5t  =  V(xi)St  +  8tSx 


dV(x ) 


dx 


X=XI 


+  0(5t5x2) 


(33) 

we  see  that  we  can  neglect  the  second  term  on  the 
RHS  because  of  diffusive  ordering  5t5x  ~  5xs  since  we 
need  to  keep  terms  only  to  order  &c2.  Therefore,  in 
the  continuum  limit  (32)  is  well  ^^roximated  by 

m- 


■>p(xi,t8)  =  (34) 

+  e~iV(-xt)St  [il>(xi+i,t0)  +  ^(mz_i,  t0)] 

-  [ip(xw,to)  +  ■ 


Now  multiplying  through  by  ezV(.xi)St  and  expanding 
the  LHS  to  order  St2  we  have  the  following  finite- 
difference  equation: 

[l  +  iV(xi)dt\ilj(xi,t8)  =  (35) 

+  ip(xl+1,t0)  +ip{xi-i,t0) 

1  —  Z 

~  [Hxl+‘2,to)  +  1p(xi-2,t0)\- 


In  the  continuum  limit,  this  finite-difference  equation 
represents  the  Schroedinger  wave  equation  with  an  ex¬ 
ternal  potential  term 


dt  {  ]  2  5t  dx 2 


— iV  (x)ip(x,  t)+0(5x4). 


(36) 

To  confirm  the  validity  of  (36)  we  perform  the  follow¬ 
ing  numerical  simulations  that  yield  results  that  can 
be  checked  against  analytical  predictions: 


It  is  possible  to  model  an  external  potential  by  apply¬ 
ing  a  local  phase  change  to  the  system  wave  function 
[13,  14] 

ip{x,t)  — >  e~'lV^x>stip{x,t).  (31) 

The  effect  of  this  phase  change  is  to  alter  the  finite 
difference  equation  (20)  as  follows 

=  -i±ViF(*‘>5V(z/,*o)  (32) 


1.  Harmonic  oscillation  of  a  displaced  Gaussian  wave 
packet  in  a  parabolic  potential. 

2.  Quantum  tunneling  through  a  potential  barrier. 

3.1  Harmonic  Oscillator 

The  first  numerical  test  presented  here  is  the  simula¬ 
tion  of  the  behavior  of  a  wave  packet  in  an  external 
parabolic  potential.  This  is  the  well-known  problem  of 
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Figure  3:  Time  evolution  of  a  Gaussian  packet  initially  dis¬ 
placed  by  a  =  321  lattice  sites  from  the  center  of  a  parabolic 
potential  well  with  K  =  10-5.  The  width  of  the  packet  is 
a  .=  14.41.  The  time  development  of  the  Gaussian  packets 
over  plotted  in  succession  where  the  x-axis  is  the  position  on 
a  L  =  256!  node  lattice  and  the  y-axis  is  the  probability 
density  \iji(x,t)  |2.  The  red  curve  is  the  parabolic  potential. 
The  h  —  1  and  m  =  1,  the  time  period  of  the  oscillation  is 
Tc  =  =  1986. 92r.  A  total  of  ten  profiles  are  over  plotted 

corresponding  to  time  t  =  0,  100r,  200r,  ...,  1000r,  which  is  ap¬ 
proximately  half  of  the  oscillation  time  period,  so  the  packet 
is  seen  to  “swing”  to  the  other  side  of  the  potential  well  while 
maintaining  a  fixed  shape  as  analytically  predicted. 

the  linear  harmonic  oscillator.  Schroedinger  analyti¬ 
cally  calculated  the  exact  time-dependent  solution  for 
the  evolution  of  a  Gaussian  packet  that  is  displaced  by 
a  distance  a  from  its  central  ground  state  in  a  parabolic 
potential  well  of  the  form  V(x )  =  ~Kx? .  The  initial 
wave  function  is 

iP(x,0)  =  ,  (37) 


Time  Step  Iteration  (t) 


Figure  4:  A  comparison  between  the  analytical  and  numerical 
predictions  of  the  location  of  an  oscillating  Gaussian  packet  in 
a  harmonic  parabolic  potential  well.  The  solid  curve  is  the  an¬ 
alytical  prediction  and  the  black  circles  are  the  numerical  data 
taken  from  the  quantum  lattice  gas  simulation  presented  in  Fig¬ 
ure  3.  In  the  simulation,  the  packet  is  initially  displaced  32  lat¬ 
tice  units  from  the  center  of  the  grid  at  lattice  node  128  for  a 
periodic  system  with  a  total  of  L  =  256!  nodes.  The  numerical 
predictions  are  in  excellent  agreement  with  the  exact  analytical 
solution. 

wave  packet  was  recorded  every  lOOr  time  steps.  This 
data  is  plotted  in  Fig  4.  The  location  of  the  peak 
oscillates  in  time  as  expected.  Overplotted  on  this 
numerical  data  is  the  exact  solution  for  the  oscillation 
a  cos  cjct+x0  and  the  agreement  between  the  analytical 
solution  and  the  numerical  data  is  excellent. 


where  a  =  (mK/h1) 4  is  the  width  of  the  packet  and 

wc  =  (K/rri)*  is  the  angular  frequency  of  the  classical 
harmonic  oscillator  [21].  The  exact  time-dependent 
solution  for  the  probability  density  is  the  following: 

\lp{x,  f)|2  =  -2Le-K(*-«cosWci)2_  (3g) 

7T2 

A  derivation  of  the  result  (38)  is  also  presented  by 
Schiff  [22], 

To  test  the  quantum  lattice  gas  algorithm  against 
(38)  we  used  a  periodic  lattice  with  L  =  256!  nodes. 
The  initial  Gaussian  packet  is  displaced  to  the  right 
of  the  center  of  the  grid  by  32  lattice  nodes  and  so 
is  initially  located  at  xa  =  160!  as  shown  in  Fig  3. 
With  h  =  1  and  m  =  1,  the  classical  time  period 
is  Tc  —  2n/ioc  =  1987t.  So  letting  the  simulation 
rim  for  1000  iterations  allows  the  packet  to  the  other 
side  of  the  potential  well  near  position  x  =  961!  as 
demonstrated  in  Fig  3. 

The  simulation  was  rum  for  a  total  of  6000  time 
steps  and  the  location  of  the  peak  of  the  Gaussian 


3.2  Scattering  Off  A  Potential  Barrier 

The  next  numerical  test  of  the  quantum  lattice  gas  is 
to  simulate  the  well-known  case  of  quantum  tunneling 
through  a  constant  potential  barrier  of  width  a.  That 
is,  V (x)  =  V0  for  0  <  x  <  a  and  V (x)  =  0  otherwise. 
The  initial  wave  function  is  a  Gaussian  packet  with 
net  momentum  to  the  right 

1>(x,  0)  =  (39) 

7T4CT4 

where  p  is  the  momentum  parameter.  We  choose  the 
mean  kinetic  energy  of  the  packet  to  be  equal  to  the 
constant  energy  level  of  the  potential  barrier  |p2  =  Va. 
In  this  case,  the  packet  tunnels  through  the  barrier 
but  the  sum  of  the  transmission  and  reflection  proba¬ 
bilities  are  less  than  one  because  there  is  a  resonance 
effect  where  the  particle  is  also  trapped  inside  of  the 
barrier.  This  effect  is  observed  in  the  numerical  simu¬ 
lation  shown  in  Figure  5. 
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4  Two  Fermionic  Particles 
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Figure  5:  A  sequence  of  snap  shots  of  the  time  evolution  of 
a  packet  that  is  incident  from  the  left  on  to  a  potential  barrier 
where  the  mean  kinetic  energy  of  the  packet  equals  the  energy 
of  the  barrier.  The  x-axis  is  the  lattice  position  and  the  y- 
axis  is  the  probability  density.  The  iteration  time  step  for  each 
frame  of  the  sequence  is  labeled  in  the  upper  left  corners.  The 
simulation  was  run  on  a  periodic  grid  of  size  L  =  4000t  for  a 
total  of  20000  time  steps.  The  width  of  the  incident  packet  was 
set  to  <j  =  .035L  =  140f  and  the  initial  momentum  parameter 
was  set  to  p  =  0.1  in  units  where  l  =  1,  r  =  1  and  m  = 
1.  The  width  of  the  barrier  was  set  to  a  =  0.064T  =  256f . 
As  expected  the  numerical  simulation  clearly  demonstrates  the 
resonance  effect  where  there  is  a  non-zero  probability  of  the 
particle  to  be  trapped  within  the  barrier  itself. 


The  efficiency  of  the  quantum  algorithm  (12)  becomes 
evident  when  it  is  used  to  simulate  the  dynamics  of 
multiple  quantum  ’particles.  The  case  of  multiple 
quantum  particles  is  still  handled  by  the  same  evolu¬ 
tion  operator  E  that  we  tested  for  the  single  particle 
case.  The  particular  sequence  and  number  of  quan¬ 
tum  gate  operations  remains  fixed,  independent  of  the 
number  of  particles  to  be  simulated.  The  only  differ¬ 
ence  is  how  the  system  wave  function  is  initialized. 

In  this  section,  for  pedagogical  reasons,  we  will  con¬ 
sider  the  case  of  simulating  two  free  quantum  particles. 
The  approach  we  use  in  this  case  can  be  directly  gen¬ 
eralized  to  the  many-particle  case. 

To  begin  with  we  write  the  Schroedinger  wave  equa¬ 
tion  for  two  free  quantum  particles 


ihdip(x,y,t)  _  ti 2  d2ip(x,y,t)  _  h2  d2ip( x,y,t) 
dt  2m  dx2  2m  dy2  ’ 

(40) 

where  x  and  y  are  the  spatial  coordinates  of  the  first 
and  second  particle,  respectively.  Since  the  wave  func¬ 
tion  is  spatially  separable  as  ip(x,  y,  t)  =  tp{x,  t)tp{y,  t), 
the  analytical  solution  to  (40)  is  obtained  in  a  similar 
fashion  to  the  one-body  case  by  computing  the  Fourier 
components  of  the  energy  basis  functions 


do 

dn 

bn 


!  fL/2 

=  j  <p(x,0)dx  (41) 

L  J-L/2 

2  fL> 2  /'  X\ 

=  —  /  ip(x,  0)  cos  \  2mv—J  dx  (42) 

2  fL /2  /  x  x 

=  -jr  J  ip(x,  0)  sin  \  2nir—j  dx.  (43) 


The  energy  eigenvalues  are  still  given  by  (27)  and  the 
time-dependent  single-particle  solution  is 


y>(x,  t) 


+ 


(44) 


which  is  basically  the  same  as  (28)  except  that  we  had 
to  add  the  sin  term  because  with  two  particles  the 
wavefunction  is  not  even,  as  is  (24)  for  example.  We 
shall  test  the  time  evolution  of  two  Gaussian  packets. 
The  initial  wave  function  in  our  test  is  the  odd  function 


tpeX«ct(x,  y,  t) 


—  ^2  [^al>al  ^)v3a2,o'2  {Vi  t) 

~  <Pai,<ri{yit')(Pa2,<r2(xyt)\’  (45) 
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where 

1  ■' 

^Pcc,(y  (*r  ,  0)  Y  re  2a~  *  (45) 

U  2  7T? 

The  subscripts  on  the  function  ipat&  denote  its  de¬ 
pendence  on  the  position  and  width  of  the  individual 
Gaussian  packet.  This  functional  dependence  is  actu¬ 
ally  contained  within  the  form  of  the  coefficients  a0, 
an ,  and  bn  that  depend  on  the  position  and  width  of 
the  Gaussian  packet  in  accordance  with  (41)  through 
(43).  Note  that  given  the  form  of  (45),  -ipexact{x,x,t)  = 

0  which  satisfies  the  Pauli  exclusion  principle. 

4.1  Numerical  Confirmation 

To  numerically  simulate  the  evolution  of  the  two- 
particle  wave  function  governed  by  (40)  using  quan¬ 
tum  algorithm  (12)  we  must  use  a  new  computational 
formulation  to  implement  our  algorithm.  The  finite- 
difference  equation  implementation  that  we  used  in 
Section  2.4,  in  the  single-particle  case,  cannot  be  di¬ 
rectly  applied  in  the  two-particle  case  to  each  particle 
individually  because  it  does  not  allow  for  the  possi¬ 
bility  when  the  particles  are  quantum  mechanically 
entangled.  In  general,  this  will  be  the  case  when  there 
is  an  interaction  between  the  particles.  Therefore,  we 
shall  use  an  implementation  that  can  handle  the  most 
general  situations  involving  correlated  particles  and 
one  that  naturally  scales  to  handle  an  arbitrarily  large 
number  of  particles  in  the  system. 

We  shall  represent  the  basic  quantum  gate  opera¬ 
tions  in  terms  of  the  fermionic  creation  and  annihila¬ 
tion  operators  in  the  number  representation,  denoted 
ah  and  aa  respectively,  and  use  this  approach  as  the 
basis  for  a  general  computational  formulation  appli¬ 
cable,  in  particular,  to  our  algorithm  and,  in  general, 
to  any  quantum  algorithm.  Acting  on  a  system  of  Q 
qubits,  a£,  and  aa  create  and  destroy  a  particle  occu¬ 
pancy  encoded  in  the  ath  qubit 

v  f  0  ,  na  =  1 

aa\ni...na...nQ)  -  j  |ni . . .  i . .. nQ)  ,  nQ  =  0 

(47) 

aa\nl...na...nQ)  -  j  0  ,  na  =  0  ' 

.(48) 

The  fermionic  creation  and  annihilation  operators  sat¬ 
isfy  the  anti-commutation  relations 

{&cn  =  (49) 

{4>4>  -  °- 


The  number  operator  na  =  a^aaa  has  eigenvalues  of 
1  or  0  in  the  number  representation  when  acting  on 
a  pure  state,  corresponding  to  the  ath  qubit  being  in 
state  |1)  or  |0)  respectively. 

The  square-root-of-swap  gate  (5)  acting  on  the  on¬ 
site  qubits  indexed  by  a  and  a  +  1  can  be  expressed 
in  terms  of  the  creation  and  annihilation  operators  as 

&a,a+l  ~  4  Tia(  1  ha+l)  Ao}aQ.a+i  ~ 

+  A*(l-na)na+i  +  l-ha-na+1,  (50) 

where  A  =  \  Also,  the  swap  operator  (9)  acting 
between  the  first  qubits  indexed  by  a  and  /?  at  neigh¬ 
boring  nodes  can  be  expressed  in  terms  of  the  creation 
and  annihilation  operators  as 

Xa,l3  =  1  -  0*0/ 3  -  ~  fla  -  ftp.  (51) 

The  quantum  gates  (50)  and  (51)  are  used  to  imple¬ 
ment  the  quantum  lattice  gas  collision  and  streaming 
operations,  respectively  [23]. 

The  basis  state  in  the  two-particle  sector  can  be 
labeled  with  the  binary  encoding  formula  |2“_1+2^-1) 
where  the  integers  a  and  /?  are  in  the  ranges  1  <  a  <  Q 
and  a  + 1  <  /?  <  Q.  The  number  of  basis  states  in  this 

case  is  the  binomial  coefficient  (  ^  )  •  ^he  system  ket 

can  then  be  expressed  in  the  two-particle  sector  as 

Q  Q 

w  =  E  E  wr'+z^1).  (52) 

a— l p—a+1 

Since  there  are  two  qubits  per  site,  we  initialize  the 
wave  function  using  (45)  as  follows: 


(l^J 


L+  1  ,/3  +  l, 

2  ’  .2  J 


_ i  <*  +  1  i  L+  1  ,  /?  +  1 1  L  +  1  ^ 

S a,/3  —  r  exact  (  L  2  -I  2  *  2  2  *  y  * 

(53) 

where  the  notation  [x\  means  the  floor  of  x  and  where 
Q  =  2L.  The  floor  operation  is  used  so  that  the  ini¬ 
tial  value  of  the  wave  function  at  each  node  is  divided 
evenly  between  each  pair  of  on-site  qubits.  This  is 
on  account  of  definition  (4)  that  allowed  us  to  have 
interfering  possibilities  for  a  single  particle  to  occupy 
a  single  position  on  the  lattice.  Moreover  in  the  two 
particle  case,  still  only  a  single  particle  can  occupy  a 
single  position  because  of  the  form  of  the  wave  func¬ 
tion  (45)  which  is  consistent  with  the  anti-commutater 
relations  (49).  However,  particle  one  can  interfere  on¬ 
site  with  itself  or  with  particle  two,  or  vice  versa  since 
the  particles  are  indistinguishable. 

At  this  point  we  have  described  how  we  implement 
the  two  quantum  gates  used  in  our  algorithm,  how 
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Position 

Figure  6:  Time  evolution  of  two  fermionic  particles  initialized 
as  Gaussian  packets  overplotted  in  succession  where  the  x-axis  is 
the  position  on  a  30-node  lattice  in  units  of  the  lattice  spacing  £ 
and  the  y-axis  is  the  probability  density  i/jfxi,  X2,  t)|2  projected 
onto  the  au-axis.  The  solid  curves  are  the  exact  analytical  solu¬ 
tion  and  the  circles  are  the  data  from  the  quantum  lattice-gas 
simulation  (the  initial  wave  function  was  normalized,  therefore 
the  area  under  each  curve  is  one).  The  initial  Gaussian  pack¬ 
ets  of  width  <r  =  3£  at  t  =  0  of  the  first  and  second  particle 
is  centered  at  x  =  I  Of  and  x  =  20t,  respectively.  The  disper¬ 
sion  of  both  packets  is  evident  by  observing  the  wave  function 
at  the  subsequent  times  t  =  7r,  21r,  28r,  35r  and  42r.  Periodic 
boundary  conditions  were  used  and  nmax  ='  40  energy  eigen- 
modes  were  used  to  generate  the  exact  solutions  at  four  times 
the  resolution  of  the  numerical  solution.  No  time  scale  factor 
was  used  and  there  is  good  agreement  between  the  analytical 
and  numerical  predictions  at  all  later  times  of  the  numerical 
simulation  as  demonstated  by  the  graphs. 


we  enumerate  the  basis  states,  and  how  we  initialize 
the  two-body  wave  function  in  this  basis.  The  only 
remaining  issue  left  to  describe  is  how  we  project  the 
two-coordinate  wave  function  ip(x,  y,  t )  on  to  a  single¬ 
coordinate  wave  function  ip{x,t)  that  can  be  plotted 
on  a  single  physical  axis.  Because  of  the  underlying 
lattice  in  our  system  this  is  straightforward  to  do  by 
summing  out  one  of  the  coordinates  as  follows: 

L-l 

ip(xh  tn)  =  ^  tp{xi,ym,  tn).  (54) 

m— 0 

If  ip{xi,ym,tn)  is  normalized  then  so  is  i/)(xi,tn )  ac¬ 
cording  to  (54).  A  comparison  of  the  time  evolution  of 
the  analytical  solution  (45)  and  the  numerical  solution 
(54)  for  a  lattice  with  30  nodes  is  shown  in  Figure  6. 
Even  with  this  small  lattice,  throughout  the  time  evo¬ 
lution  of  the  model  run  the  numerical  predictions  are 
in  good  agreement  with  the  predictions  of  the  exact 
solution. 


5  Conclusion 

We  have  presented  a  quantum  algorithm  that  is  an 
efficient  and  accurate  scheme  for  simulating  the  time- 
dependent  evolution  of  a  system  of  quantum  particles 
governed  by  the  non-relativistic  Schroedinger  wave 
equation.  The  scheme  uses  a  quantum  lattice  gas 
system  of  particles  colliding  and  hopping  on  a  lat¬ 
tice.  The  algorithm  is  efficient  in  the  sense  that  the 
computational  effort  needed  to  simulate  an  arbitrar¬ 
ily  large  number  of  particles  (within  the  constraint  of 
the  grid  resolution)  exactly  equals  the  computational 
work  needed  to  simulate  a  single  particle,  given  that 
the  algorithm  is  executed  on  a  quantum  computer  that 
remains  phase-coherent  throughout  the  entire  coarse 
of  the  simulation.  However,  a  limitation  does  exist 
on  state  preparation  [i.e.  initialize  the  quantum  com¬ 
puter’s  memory)  and  we  have  not  argued  here  that 
it  is  possible  to  initialize  the  many-body  wave  func¬ 
tion  in  an  efficient  way.  For  that  matter,  have  we  also 
have  not  argued  that  it  is  possible  to  measure  the  fi¬ 
nal  state  (reading  the  quantum  computer’s  memory)  of 
the  computed  wave  function  in  an  efficient  way.  Never¬ 
theless,  the  quantum  algorithm  presented  here,  which 
is  a  way  of  representing  a  discretized  Feynman  path 
integral,  has  the  useful  feature  that  it  is  explicit  in 
time  where  the  value  of  the  wave  function  at  location 
x  and  time  t+r  depends  only  on  the  previous  values  of 
the  wave  function  at  time  t  in  the  immediate  vicinity 
of  x.  Since  the  algorithm  is  unitary  and  fourth-order 
accurate  in  space,  it  is  useful  even  for  implementation 
on  a  classical  computer. 

We  have  carried  out  a  variety  of  numerical  tests 
proving  that  the  quantum  algorithm  indeed  allows  us 
to  faithfully  reproduce  the  correct  dynamical  behavior 
of  a  continuous  and  differentiable  wave  function  in  the 
presence  of  an  external  potential.  However,  the  total 
probability  for  finding  the  quantum  particle  in  the  sys¬ 
tem  is  not  exactly  conserved  in  this  quantum  lattice 
gas  model.  One  must  approach  the  continuum  limit 
to  achieve  a  high  degree  of  probability  conservation. 
Since  we  have  demonstrated  that  the  quantum  lattice 
gas  model  is  fourth-order  convergent  in  space,  it  is  al¬ 
ways  possible  to  choose  a  grid  resolution  that  achieves 
the  necessary  fixed  numerical  accuracy  required  by  any 
application. 

We  have  also  described  and  carried  out  the  numer¬ 
ical  simulation  of  two  fermionic  particles,  which  are 
non-interacting  except  for  a  quantum  mechanical  ex¬ 
change  force  arising  from' the  anti-commutation  rela¬ 
tions.  The  numerical  formalism  used  to  implement 
the  quantum  lattice  gas  algorithm  represents  the  basic 
quantum  gates  in  terms  of  quadratic  products  of  cre- 
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ation  and  annihilation  operators  in  a  second  quantized 
representation.  This  formalism  straightforwardly  ,han- . 
dies  an  arbitrarily  large  number  of  particles.  The  sim¬ 
ulation  is  carried  out  in  the  many-body  sector  (he. 
a  Fock  space  where  the  number  of  particles  is  fixed) 
where  all  the  basis  states  are  enumerated  by  a  sim¬ 
ple  binary  encoding  formula..  In  general,  the  wave 
function  is  initialized  using  a  Slater  determinant  so 
that  the  Pauli  exclusion  principle  is  satisfied  and  the 
wave  function  is  odd.  In  all  the  test  cases  (single  free 
particle,  single  particle  in  an  external  potential,  and 
two  free  fermions)  the  numerical  predictions  agreed 
extremely  well  the  analytical  predictions  and  exact  so¬ 
lutions. 

It  is  important  to  realize  that  it  is  possible  to  push 
the  quantum  lattice  gas  model  into  regimes  where  the 
numerically  predicted  results  are  absolutely  wrong. 
This  occurs  when  the  local  configurations  are  far  from 
local  equilibrium.  Remember  that  local  equilibrium 
exists  on  a  lattice  node  when  the  qubits  at  that  node 
have  identical  phase.  Therefore,  large  gradients  in 
the  macroscopic  profile  of  the  modeled  wave  function 
may  cause  a  large  phase  difference  between  the  on-site 
qubits  eventually  resulting  in  anomalous  large-scale 
behavior  in  the  model.  A  good  example  of  this  occurs 
when  the  momentum  of  a  moving  wave  packet  is  too 
high.  In  this  case,  since  the  real  and  imaginary  parts 
of  the  wave  function  are  sinusoidal,  if  the  momentum 
is  so  high  that  the  wavelength  of  the  traveling  wave 
is  on  the  order  of  the  lattice  cell  size  (A  =  %/p  ~  £), 
then  after  a  few  times  step  iterations  of  the  algorithm, 
large  phase  differences  in  the  on-site  qubits  occur.  In 
this  case,  the  local  on-site  configuration  will  not  re¬ 
lax  toward  the  correct  local  equilibrium.  Fortunately, 
the  norm  of  the  modeled  wave  function  will  deviate 
from  unity  in  these  types  of  cases.  Therefore,  it  is 
straightforward  to  test  if  the  model  predictions  are 
non-physical  by  periodically  checking  the  norm  of  the 
wave  function. 

Interaction  potentials  between  particles  can  also  be 
modeled  using  this  quantum  lattice  gas  method  [14]. 
We  still  need  to  perform  simulations  of  a  many-body 
system  with  an  interaction  potential.  Also,  numerical 
tests  in  two  and  three  dimensions  should  be  conducted. 
In  this  case,  it  would  be  straightforward  to  test  the  ad¬ 
dition  of  an  external  vector  potential  using  the  results 
analytically  predicted  by  Polley  [18].  The  simulation 
of  the  dynamical  behavior  of  the  positronium  atom 
would  be  a  reasonable  next  step  for  the  application  of 
the  quantum  lattice  gas  method  to  quantum  mechan¬ 
ical  computational  physics. 
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A  Appendix 

The  full  finite-difference  equation  for  the  quantum 
lattice-gas  model  presented  in  this  paper  is  very  long. 

To  simplify  this  expression,  we  introduce  a  local  neigh¬ 
borhood  vector  with  the  following  18  components 

(55) 

<Po(xi+l,  t),  ip1{xi+1,t),  <po(xi-i,t), 

(po (x{+2 ,  t) ,  <pi  (xi+2 ,t),<p0{xi-2,t),  ipi  (xi-2,t), 
<fi0(xi+3,t),  (fii  (xt+ 3,  t ),  (p0(xi-3,  t),  ip i  (xi-3,  t ), 
(p0(xi+i,  t),  (pi(xi+4 ,  t),  <p0(xi-4,t),  <fil(xi-4,  t )). 
We  define  the  following  two  coefficient  vectors 
a  =  (— 3,  —  3i,  3,  i,  —  3,  —  5i,  1,  i,  3, 3i,  —  1,  i,  1, 3i,  0, 0,  —  1,  —  i) 
(3  =  (— 3i,  —  3,  —  5i,  —  3,i,  3, 3i,3,i,  1, 3i,  l,t,  —  1,  —  i,  —  1, 0, 0). 

(56) 

The  microscopic  evolution  equation  (13),  explicitly 
written  out,  has  the  following  protocol  of  operations 

Wtie)>  =  (s^CS2CSfCS2C)  (slC^CSlCSiC)  iV’(to))- 

(57) 

The  corresponding  full  finite-difference  equation  can 
be  specified  by  the  following  dot  product  of  these  vec¬ 
tors 

Po(xi,tie)  =  -  (58) 

Vi  K<„)  =  (59) 

Note  that  if  to  is  the  initial  time,  then  the  interval 
r  =  tie  is  defined  the  update  time  step.  The  finite- 
difference  equation  for  ip  =  ipo  +  ipi  is 

«*»(,.)  (60) 

and  it  has  a  high  degree  of  numerical  accuracy  as  in¬ 
dicated  in  Figure  2. 
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